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Abstract 
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Although  a  number  of  the  major  objectives  that  were  established  at 

/ 

^the  outset  of  this  project  were  not  met, /a  number  of  milestones  were 
realized r^The  digital  implementation  of  a  negating  phase  shift  that 
operates  perfectly  under  ideal  conditions  was  a  major  accomplishment. 


The  establishment  of  a  zero  level  of 


:ions  wj 
d 


was  also  significant.  The 


incorporation  of  the  exponential  smoothing  technique  to  minimize  the 
effect  of  measurement  noise  was  important  since  it  uncovered  a  possible 
connection  between  the  size  of  the  target  image  and  its  performance 
throughout  the  pattern  recognition  process.  However,  the  major  obstacle 
that  surface  during  the  execution  of  this  project  was  a  filter  diver¬ 
gence  problem.  It  has  been  proposed  that  this  problem  can  be  solved  by 
implementing  the  Fourier  transform  derivative  property  instead  of  the 
forward-backward  difference  method  to  compute  the  spatial  derivative  of 
the  non-linear  h  function. 
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L  Introduction 


Background 

The  application  of  laser  technology  to  everyday  life  is  growing  in 
importance  as  each  year  passes.  In  areas  from  medicine  to  industry  to 

j 

military  applications,  laser  technology,  although  in  its  embryonic  stage, 
has  gained  a  foothold  and  is  destined  to  have  a  major  impact  on  the 
future.  For  instance,  in  the  industrial  area,  laser  drilling  has  signi¬ 
ficantly  improved  the  machining  qualities  resulting  in  smoother  cuts, 
reduced  tool  force,  and  increased  speed  and  accuracy  (1:225-31).  In  the 
medical  area,  the  application  of  lasers  in  opthalmology  is  promising. 

For  example,  the  use  of  lasers  to  repair  retinal  tears  and  holes  has  been 
quite  successful  (2:360-7).  For  military  applications,  the  ability  to 
deposit  large  amounts  of  laser  energy  onto  targets  is  a  major  research 
effort.  One  in  a  number  of  major  obstacles  before  realizing  this  appli¬ 
cation  deals  with  the  precision  tracking  of  a  target  and  the  subsequent 
pointing  of  the  laser  beam.  This  thesis  is  one  in  a  series  of  research 
efforts  devoted  to  solving  this  problem. 

Traditionally,  correlation  algorithms  have  been  employed  to  provide 
pointing  and  tracking  information  to  a  system.  This  correlation  tracker 
stores  a  set  of  predetermined  or  previous  real-time  data  target  images 
in  memory,  compares  these  images  with  the  images  from  its  sensors.  In 
the  process  of  performing  the  mathematical  calculations  to  characterize 
the  differences  between  the  stored  and  actual  images,  the  appropriate 
commands  are  sent  to  the  tracker  to  minimize  the  offset  between  the  two 
(3:30-8).  However,  two  major  disadvantages  of  the  correlation  algorithm 
are  its  suspectibility  to  noise  and  absence  of  sensitivity  to  target 
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dynamics  (3:31).  To  combat  those  drawbacks,  the  use  of  an  extended 
Kalman  Filter  algorithm  in  place  of  the  correlation  tracker  ..is  being 
explored. 

The  Kalman  Filter  is  a  computer  algorithm  which  processes  noise- 
corrupted  measurements  and  provides  a  reasonably  accurate  estimate  of 
the  state  variables  of  interest  (4:3).  The  actual  mathematical  details 
of  the  Ka '.man  Filter  are  contained  in  Chapter  2  under  the  subtitle 
Extended  Kalman  Filter. 

As  mentioned  earlier,  this  project  is  one  in  a  series  of  research 
efforts.  To  be  more  precise,  this  project  is  the  follow-on  to  two  other 
projects.  The  first  of  these,  entitled  "An  Extended  Kalman  Filter  For 
Use  in  a  Shared  Aperture  Medium  Range  Tracker"  by  Daniel  E.  Mercier, 
dealt  with  a  very  benign  distant  point  target  which  could  be  analytically 
expressed  as  a  two-dimensional  gaussian  intensity  profile  of  circular 
contours  (3:6-7).  This  intensity  profile  and  all  other  profiles  devel¬ 
oped  for  this  research  effort  are  assumed  to  be  scanned  by  a  Forward- 
Looking  Infrared  (FLIR)  sensor.  This  sensor  horizontally  and  vertically 
scans  the  system  field  of  view  (FOV)  to  provide  an  8  x  8  array  of  dis¬ 
crete  values.  These  discrete  values  represent  the  average  intensity  of 
that  portion  of  the  image  which  lay  across  the  appropriate  detector  (3: 
4-5). 

Using  some  of  the  results  of  Mercier' s  Thesis,  the  second  project, 
entitled  "An  Adaptive  Distributed  Measurement  Extended  Kalman  Filter  For 
a  Short  Range  Tracker"  by  Robert  L.  Jensen  and  Douglas  A.  Harnly,  dealt 
with  more  dynamic  targets  at  closer  ranges  but  still  assumed  that  the 
target  intensity  pattern  could  be  expressed  analytically.  However, 

Jensen  and  Harnly's  thesis  did  make  provisions  to  adaptively  change  the 
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shape  of  unimodal  intensity  patterns  (5:77).  Still,  for  actual  hardware 
implementation,  a  priori  knowledge  of  the  analytic  form  of  the  intensity 
often  cannot  he  assumed.  Instead,  on-line  numerical  techniques  will 
have  to  be  exploited  to  provide  the  necessary  information  to  the  Kalman 
Filter.  The  development  of  these  techniques  and  their  subsequent  inter¬ 
face  with  the  Kalman  Filter  is  the  basis  for  this  project. 

Problem  Overview 

The  numerical  techniques  mentioned  toward  the  end  of  the  previous 
section  include  the  Fast  Fourier  Transform  (FFT),  the  shift  theorem  of 
Fourier  Transform,  the  exponential  smoothing  technique,  arid  the  forward- 
backward  difference  method.  These  techniques  are  discussed  in  more 
detail  in  Chapter  2  and  Appendix  A.  However,  all  involve  the  manipula¬ 
tion  of  intensity  measurements  provided  by  the  FL1R  sensor.  The  end 
result  of  these  computations  is  to  provide  information  to  the  Kalman 
Filter  about  the  intensity  pattern  shape.  This  information  takes  the 
form  of  certain  components  of  the  measurement  update  equation.  A  de¬ 
tailed  description  of  the  measurement  update  equation  is  contained  in 
Chapter  2  under  the  subtitle  Measurement  Update.  However,  the  general 
form  of  this  equation  is 

x(t+)  =  x(t~)  +  K(t.)  (z(t.)  -  h(x(t-)t.)  (1) 

where  £(t"£)  ~  state  estimate  vector  after  incorporation  of  measurement 

x(t~)  *  propagated  state  estimate  vector  before  incorporation  of 
measurement 

K(t^)  =  Kalman  Filter  Cain 

2_(t^)  -  actual  measurement  vector 
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h(x(tT)t,))  =  non-linear  function  of  intensity  measurement  at 

—  i  1 

time  t^,  as  a  function  of  the  true  state  osLimatc 

In  devising  this  project,  it  was  decided  that  under  ideal  conditions,  the 
Kalman  Filter  would  provide  state  estimates  which  would  center  the  target 
images  from  one  sample  period  to  another.  This  desirability  of  producing 
centered  images  is  motivated  by  the  simplicity  it  produces  in  the  Kalman 
Filter  equations.  If  the  center  images  correspond  to  state  estimates 
equal  to  zero,  the  resulting  measurement  equations  become 

i(tT)  =  0  (2) 

x(t+)  =  0  +  K(t.)  (z(t.)  -  h(01t.))  (3) 

where 

h^O^tJ  =  centered  non-linear  h  function 

The  generation  of  involves  the  use  of  the  FFT,  shifting  theorem, 

and  exponential  smoothing  techniques  mentioned  earlier.  The  details  of 
how  these  techniques  are  utilized  will  be  discussed  later  however,  buried 
within  the  Kalman  Filter  gain  (K(t^))  is  the  spatial  derivative  of  the 
centered  non-linear  h  function  (see  Chapter  2,  Equation  56).  This 
spatial  derivative  is  generated  using  the  forward-backward  difference 
method  discussed  in  Appendix  A. 

The  flow  of  the  data  processing  scheme  is  shown  in  Figure  1.  In 
essence,  there  are  two  parallel  data  processing  paths  for  the  intensity 
measurements.  The  first  path  involves  taking  the  8x8  array  of  inten¬ 
sity  measurements  and  arranging  it  by  rows  into  a  64  x  1  measurement 
vector.  This  vector  is  then  provided,  as  a  measurement  z_(t^),  to  the 
extended  Kalman  Filter  which  in  turn  provides  state  estimates  that  are 
used  in  the  second  path  to  provide  centered  measurement  functions.  This 
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Figure  1.  Data  Process 


Its  purpose  is  to 


second  path  represents  the  main  thrust  of  this  thesis, 
provide  the  centered  non-linear  and  linear  measurement  functions  mentioned 
earlier.  To  accomplish  this,  the  shifting  theorem  of  Fourier  Transforms 
i  exploited.  In  essence,  the  translation  of  an  intensity  pattern  in  Che 
space  domain  can  be  negated  by  multiplying  its  Fourier  Transform  by  the 
complex  conjugate  of  the  resulting  linear  phase  shift  (see  Chapter  2, 
Fourier  Transform) .  The  source  of  the  image  translation  about  the  FL1R 
array  can  be  traced  to  two  effects:  the  actual  target  dynamics  and  atmos¬ 
pheric  jitter.  An  estimate  of  these  effects  are  available  from  the  Kalman 
Filter  measurement  update  equations  (see  Chapter  2,  Extended  Kalman  Fil¬ 
ter).  In  tuvn,  these  best  estimates  arc  used  in  the  argument  of  the 
complex  conjugate  of  the  linear  phase  shiLL  to  provide  centered  measure¬ 
ment  functions.  In  addition,  before  the  inverse  Fourier  Transform  is 
taken,  this  centered  pattern  is  averaged,  using  the  exponential  smoothing 
technique,  with  previous  centered  Fourier  Transformed  patterns  to  mini¬ 
mize  the  effect  of  measurement  noise  (see  Chapter  2,  Averaging).  Once 
the  inverse  Fourier  Transform  is  taken,  the  spatial  derivative  is  then 
taken  using  the  forward-backward  difference  approximation  discussed  in 
Appendix  A  to  provide  the  centered  linearized  H  function  used  in  the 
calculation  of  the  Kalman  Filter  Cain. 

The  sequential  processing  along  the  second  path  shown  in  Figure  1 
together  with  its  connection  to  the  first  path  is  an  LmporLant  fact.  A 
copy  of  the  8x8  array  of  data  is  first  padded  with  zeroes  to  alleviate 
problems  in  using  the  l'FT  (see  Chapter  2,  Fourier  Transform).  Next,  the 
Fourier  Transform  of  the  two-dimensional  array  of  data  is  taken,  before 
the  negating  phase  shift  is  applied  to  this  Fourier  Transform,  the  pro¬ 
cessing  of  the  measurement  vector  along  the  first  path  has  to  be 
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completed,  since  the  result  ing  state  estimates  from  tlie  measurement  up¬ 
date  eciuation  arc  used  in  the  linear  phase  shift  La  the  second  data  pro¬ 
cessing  path.  Alter  applying  the  negating  linear  phase  shift,  this 
centered  l-'ourier  Transform  is  averaged  with  past  centered  Fourier  Trans¬ 
forms  to  minimize  the  effect  of  noise.  Next,  the  inverse  Fourier  Trans¬ 
form  is  taken  which,  theoretically ,  results  in  a  centered  pattern  with 
noise  effects  substantially  reduced.  The  zeroes  that  were  initially 
added  to  the  input  array  are  stripped  away  with  the  remaining  8x8  data 
array  representing  the  non-linear  h  function.  In  addition,  using  the 
numerical  approximation  discussed  in  Appendix  A,  the  linearized  H  func¬ 
tion  is  generated  from  this  non-linear  h  function  and  both  functions  are 
used  by  the  Kalman  Filter  for  processing  of  the  next  measurement  that 
becomes  available. 

An  additional  process  that  occurs  along  the  first  path  is  the  pro¬ 
pagation  of  the  state  estimate  vector  to  the  next  sample  time.  This 
process  provides  the  best  prediction  of  where  the  target  will  be  lc"?.ted 
just  before  the  next  measurement  update  is  taken.  Therefore,  this  in¬ 
formation  could  be  fed  to  a  controller  so  as  to  minimize  the  perturba¬ 
tions  of  the  image  about  the  center  of  the  FOV. 

Flan  of  Attack 

The  plan  of  attack  presented  here  provides  a  general  flow  of  what 
will  later  be  examined  in  greater  detail  in  the  performance  analysis 
section.  The  verification  of  a  pattern  recognition  algorithm  along  with 
its  interface  with  the  Extended  Kalman  Filter  represents  the  basic  pre¬ 
mise  upon  which  this  thesis  project  was  developed. 

The  verification  of  the  pattern  recognition  algorithm  is  composed 
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of  several  different  parts.  These  parts  include  verifying  the  TFT 
algorithm,  verifying  the  ability  to  negate  the  translational  effects  in 
the  space  domain  given  perfect  phase  shift  information  to  apply  in  the 
spatial  frequency  domain,  and  verifying  the  exponential  smoothing  tech¬ 
nique  as  a  means  to  minimize  the  effect  of  noise  corruption  on  the 
measurement  information. 

The  Extended  Kalman  Filter  section  involves  measuring  the  impact  of 
the  non-linear  and  linearized  h  functions  developed  in  the  pattern  rec¬ 
ognition  section  on  the  performance  of  the  Kalman  Filter  as  compared  to 
being  given  correct  h  functions.  Also  implicit  in  this  verification  is 
another  comparison  involving  the  use  of  state  estimates  from  the  Kalman 
Filter,  instead  of  providing  artificial  knowledge  of  Lhe  true  pattern 
offset,  to  provide  the  shift  information  needed  in  the  pattern  recogni¬ 
tion  section  to  provide  centered  patterns  used  to  generate  the  h  func¬ 
tions  mentioned  earlier  (see  Problem  Overview) . 


8 


II  Models  ;md  Dnlu  Processing  Fundamentals 


Karhunen-Locve  Transformation 

In  the  area  o£  pattern  recognition,  it  is  highly  desirable  to  deter¬ 
mine  the  optimal  set  of  eigenfunctions  and  their  corresponding  eigenvalues 
to  represent  two-dimensional  intensity  patterns.  This  representation  can 
be  determined  with  the  orthogonal  Karhunen-Lo'evc  transformation 

K(x,oi;y,B)  4> .  Cot, 3)  dadtT  =  \.  (a, 6)  (4) 

-y  h  x  xi 

where  R(x,a;y,3)  =  spatial  autocorrelation  kernel 
=  eigenfunctions 

=  eigenvalues 

Based  on  properties  specified  through  the  spatial  autocorrelation 
kernel,  the  entire  information  of  an  image  is  preserved  by  the  given  set 
of  eigenvalues  and  eigenfunctions  (6:6).  Furthermore,  since  the  low 
order  terms  normally  provide  the  maximum  sensitivity  to  target  motion, 
the  higher  order  terms  of  this  spatial  model  decomposition  can  be 
neglected  with  little  degradation  in  the  quality  of  the  resultant  image 
(6:6). 

As  mentioned  earlier,  the  key  to  finding  the  optimal  set  of  func¬ 
tions  is  to  know  the  spatial  autocorrelation  kernel.  For  intensity 
patterns  represented  by  a  two-dimensional  discrete  array  of  values,  the 

correlation  kerneL  is  represented  as  a  correlation  matrix  of  dimension 

2  2 

N  x  N  with  N  being  tiie  dimension  of  the  square  input  matrix  (6:6). 

The  generation  of  tills  correlation  matrix  is  a  major  drawback  in  using 
the  Karhunen-Loeve  transformation.  However,  if  this  matrix  can  be 
assumed  or  determined,  the  transformation  matrix  (A)  which  diagonalizes 
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the  correlation  matrix  (K)  can  be  evaluated  by: 


where  A^'s  represents  the  eigenvalues  in  decreasing  order  of  magnitude 
(X. >A  >  ...  >X,,2 )  (6:6).  Since  the  transformation  matrix  (A)  is  devel- 
oped  from  orthogonal  eigenvectors  of  the  correlation  matrix,  the  trans¬ 
formation  matrix  is  orthogonal.  Thus,  the  forward  and  reverse  Karhunen- 
Loeve  transformations  that  will  preserve  the  quality  of  the  image  are 
forward  direction:  (!•’)  =  (f)  (A) 

T 

reverse  direction:  (f)  =  (F)  (A) 
where  (f)  =  image  vector 

(F)  »  resultant  image  vector  in  tlie  transform  domain 
(A)  =  transforms tion  matrix  (6:7) 

Once  again,  if  the  entire  group  of  eigenvalues  and  eigenvectors  were 

retained,  the  total  quality  of  the  image  would  be  preserved.  But,  if 

2 

only  m  of  the  first  N  eigenvalues  were  retained,  the  resultant  mean 
square  error  (MSE)  between  the  reconstructed  image  and  the  initial  image 
would  be  the  sum  of  the  eigenvalues  not  included  in  the  transformation. 

N2 

MSE  -  X  X.  (6:7) 
i=m+l 

As  discussed  in  the  article  "image  Processing  by  Computer"  by  Guy 
Hanuise,  from  which  the  Karhunen-Loove  transformation  argument  has  been 
developed,  this  transformation  docs  possess  some  major  disadvantages. 

The  most  significant  drawbacks  center  around  the  generation  of  the 
correlation  matrix.  As  the  size  of  the  square  image  matrix  (N  x  N) 
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2  2  2 

grows,  the  correlation  matrix  grows  as  N  (i.e.  N  x  N  ).  Thus,  with 
large  image  quantization  levels,  serious  data  processing  problems  arise 
(6:7).  Along  the  same  lines,  the  exact  calculation  of  the  correlation 
matrix  is  very  difficult  to  perform.  As  a  result,  more  common  orthogonal 
transformations  such  as  the  Fourier  transform  are  used  in  signal  proces¬ 
sing.  Implicit  in  the  use  of  the  Fourier  transform  are  the  assumptions 
of  spatial  stationarity  of  the  autocorrelation  kernel  and  a  space  domain 
infinite  in  extent.  Referring  back  to  equation  A,  in  applying  the  assump 
tion  of  spatial  stationary,  and  a  domain  of  infinite  extent,  the  Karhunen 
Lo'eve  transformation  equation  becomes 

CO  00 

/  /  R(x-«,y-3)v^(a,  3)dad3  =  (6) 

— oo  —to 


Upon  close  examination  of  equation  0,  the  integral  equation  would  be 
recognized  as  a  two-dimensional  convolution  of  two  functions  (7:10). 

The  unique  feature  of  the  convolution  theorem  that  makes  it  very  appeal¬ 
ing  is  that  in  the  other  transform  domain,  the  two  functions  are  multi¬ 
plied  together.  Therefore,  the  Fourier  transform  of  equation  6  can  be 
written  as 

S(w  ,w  )  4>.(jw  ,iw  )  =  A.0(iw  ,w  )  (7) 

v  x’  y  i  J  x,J  y  i  J  x  y 


where  S(w  ,w  )  =  Fourier  transform  of  R(x,y)  (power  density  spectrum) 
x  y 

( iw  ,  jw  )  =  Fourier  transform  of  <»>.(u,P) 
i  x  J  y  l 

For  tlie  equality  in  equation  7  Lo  hold  either  one  of  two  conditions  must 


be  met:  either  S (w  ,w  )  -  A.,  which  is  an  impossibility  since  the  A.'s 
x  y  i  i 

are  scalar  multiples  and  S(w  ,w  )  is  in  a  functional  form,  or  more 

x  y 

realistically,  $.(jw  , jw  )  is  an  impulse  function.  The  choice  of  the 
i  x  y 

impulse  function  would  be  appropriate  since  it  can  be  set  to  sample  the 
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power  ileus  it.  V  spectrum  at  Lho  pa  r  I  i  cu  1 .1  r  Value  ul  w  1 where  the  function 

x 

equals  the  value  of  the  scalar  multiple. 


<!\  (jw  ,  jw  ) 
i  J  x,J  y 


b(w  -w  . 

X  XI 


; w  -w  . ) 

y  y  r 


(8) 


Thus,  by  taking  the  inverse  Fourier  transform  of  this  impulse  function, 
the  resulting  eigenfunctions  for  the  space  domain  are: 

$.(x,y)  =  T  ^(6(w  -w  .  ;w  -w  .))  (9) 

i  x  xi  y  yi 

*  *.(x,y)  =  exp(j(wx./.  -i-w^.X))  (8:237-42)  (10) 


In  reality,  the  space  domain  is  limited  by  the  system  FOV  and  the  random 
process  describing  tile  image  intensity  may  not  be  truly  spatially  sta¬ 
tionary.  However,  if  the  system  FOV  is  relatively  large  and  the  random 
process  is  quasi-stationary ,  i>.  could  be  hour  istically  argued  that  the 
resulting  eigenfunctions  still  asymptotically  approach  complex  exponen¬ 
tials.  This  result  would  provide  reasonable  motivation  to  use  complex 
exponentials  (Fourier  transform)  as  the  transformation  function  on  the 
images  in  question. 

In  closing  this  section,  it  should  be  reiterated  that  the  Karhunen- 
Loeve  transformation  is  difficult  to  perform  in  its  exact  form.  However, 
under  the  assumptions  of  spatial  stationarity  and  a  space  domain  large  in 
extent,  the  Karhunen-Loeve  equation  provides  adequate  motivation  for 
using  familiar  transformations  involving  complex  exponentials  such  as 
Four  ie  r. 


*  Tl: is  same  argument  has  been  developed  for  the  one-dimensional  case  in 
Chapter  8  of  Introduction  to  Statistical  Pattern  liecogn  i  t  ion  by  Kcinosuk.e 
Fu kuna go . 
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Fourier  Trans  to rm 


The  Fourier  transform  is  a  familiar  transformation  to  the  electrical 
engineer.  In  the  one  dimensional  case,  it  is  a  transform  quite  often 
used  to  relate  occurrences  in  the  time  domain  to  those  in  the  frequency 
domain.  However,  with  the  ability  of  lenses  to  perform  Fourier  trans¬ 
forming  instantaneously,  the  field  of  Fourier  optics  has  provided  moti¬ 
vation  for  extending  the  concept  of  the  Fourier  transform  and  its  prop¬ 
erties  into  two-dimensional  space.  As  a  result  of  this  extension,  the 
two-dimensional  Fourier  transform  becomes 

OO 

u(f  ,f  )  =  /  f  g(x,y)  exp  —  j 2rr ( f  x  +  f  y)  dxdy  (7:5)  (11) 

x  y  _  jo  x  y 


..here 

G(f  , f  )  =  Fourier  or  Frequency  Spectrum 
x  y 

g(x,y)  =  Function  in  the  Space  Domain 

f  , f  =  Soatiai  Frequencies 
x  y 

x,y  =  Spatial  Variables 


In  comparing  the  one-dimensional  and  two-dimensional  Fourier  trans¬ 
forms,  similarities  should  be  recognizable;  the  use  of  the  complex 
exponential  as  the  eigenfunction,  and  the  generation  of  spatial  frequen¬ 


cies  f  ,f  to  correspond  with  the  spatial  coordinate  x  and  y.  A  parti- 
x  y 


cular  property  that  has  been  extended  to  the  two-dimensional  case  which 
is  a  vital  part  of  this  thesis  is  the  shift  theorem.  The  shift  theorem 
states  that  a  translation  of  an  image  in  the  space  domain  results  in  a 
linear  phase  shift  in  the  spatial  frequency  domain: 


?{g(x-a,y-b) }  =  G(f  ,f  )  exp(-j2ir(f  a  +  f  b))  (7:9)  (12) 

x  y  x  y 


where 


a  =  o  I  f  so  L  of  Lite  spaLlnl  I' line  lion  ailing  llir  x  direction  from  a 
centered  position 

b  =  offset  of  the  spatial  function  along  the  y  direction  from  a 
centered  position 

To  negate  the  translational  effects  in  the  space  domain,  the  Fourier 
transform  of  the  translated  image  is  multiplied  by  the  complex  conjugate 
cf  tl-.e  linear  phase  shift; 

g(x,y)  =  T  1{G(f  ,f  )  exp (-j 2u (f  a  +  f  b)  )cxp  ( j  2it  (f  a  +  f  U)  )  }  (13) 

x  y  x  y  x  y 

Up  to  now,  the  continuous  case  of  the  Fourier  transform  has  been 
discussed,  but  to  utilize  the  Fourier  transform  and  its  properties  for 
computer  simulation,  the  discrete  form  of  the  Fourier  equations  must  be 
used.  Such  a  development  is  contained  in  Digital  Signal  Processing  by 
Alan  V.  Gppenheim  and  Ronald  V.'.  Schafer.  The  discrete  Fourier  transform 
and  the  discrete  version  of  the  shift  theorem  are  as  follows: 

>1-1  N-l 


'(x(m ,n) )  =  Z  L  x(kAK'xp(  exp  (  (14) 

ll"  k=0  v-0 


Discrete  Fourier  Transform  (9:113) 


M-l  N-l 

y  y 

k=o  a=o 


T(x(m-mn,n-n^))  =  ~  ^  >-  X(’x,^-)cxp(  •L~--m)  exp  (  J-~-— ) 


-jil'Jdl'! 
N 


.-jZ’Tkmo.  ^  —  j  2 <-  i-nos 

x  exp (  - )  exp  (  N-  -  -) 


(15) 


0  <  m  M  0<n  <N 
—  o  —  —  o  — 


Discrete  Shift  Theorem  (9:110) 


where  m.n  =  sequence  numbers  in  the  space  domain 

k ,  C-  =  sequence  numbers  in  spatial  frequency  domain 
M.N  =  spatial  area  covered  by  the  sequence  numbers 
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ro  ,n  =  translation  of  discrete  pattern  about  its  centered 
oo 

locat ion. 

As  the  definition  of  the  variables  may  imply,  the  use  of  the  discrete 
Fourier  transform  implies  certain  knowledge  is  known.  In  essence,  the 
discrete  Fourier  transform  views  a  finite  area  as  being  repeated  in¬ 
definite  in  both  the  x  and  y  directions.  Therefore,  the  establishment 
of  the  two-dimensional  periodicity  is  imperative.  The  discretzation  of 
the  finite  area  is  reflected  in  the  sequence  numbers.  The  relationship 
between  the  sequence  numbers  and  distance  in  the  original  space  domain 
is  shown  in  Figure  2.  Notice  that  each  pair  of  sequence  numbers  repre¬ 
sent  a  smaller  area  within  the  previously  mentioned  finite  area.  For 
application  to  this  thesis,  the  smaller  area  is  represented  by  20  urad  x 
20  urad  with  the  total  finite  area  being  480  urad  x  480  urad.  This 
results  in  the  sequence  numbers  (m,n)  varying  from  1  to  24.  Thus  the 
sequence  lengths  (M,N)  equal  24. 

In  the  computer  software  for  this  project,  a  more  efficient  version 
of  the  discrete  Fourier  transform  known  as  the  Fast  Fourier  Transform 
(FFT)  is  used.  The  KFT  and  the  Inverse  Fast  Fourier  Transform  (1FFT) 
are  performed  using  a  subroutine  called  F0UKT.  F0URT  is  a  multi¬ 
dimensional  FFT  routine  which  uses  the  Cooley-Tukey  method  of  calculation 
(10:76-9).  A  unique  feature  of  the  F0URT  subroutine  is  the  arrangement 
of  the  data  array  in  Liu-  spatial  frequency  domain.  As  a  result  of  the 
FOUR!  software,  the  F0URT  format  locates  the  D.C.  term  (zero  frequency) 
in  the  upper  left  hand  corner  and  the  i.armonics  are  msaligncd  as  shown 
in  Figure  3.  In  order  to  perform  signal  processing  in  the  spatial  fre¬ 
quency  domain,  the  data  array  quadrants  arc  switched,  2  with  4  and  1  with 
.3  (11:17).  As  shown  in  Figure  4,  for  an  M  by  N  dimension  array,  where  M 
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Figure  2A.  Finite  Area  A  x  B 
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and  N  have  to  be  even  numbers,  this  rearranged  format  results  in  the 

M  N 

D.C.  term  located  at  (—  +1,  -  -  +  1)  with  the  first  concentric  window 
surrounding  the  D.C.  term  being  the  first  harmonic  and  corresponding 
higher  harmonics  matching  with  its  appropriate  window  (11:18).  Once 
the  signal  processing  has  taken  place,  the  data  array  lias  to  be  re¬ 
arranged  back  into  the  FOURT  format  before  the  IFFT  is  taken  (11:19-21). 

Another  requirement  for  FFT  processing  is  to  pad  the  input  data 
array  with  zeroes  to  reduce  edge  effects,  aliasing,  and  leakage  condi¬ 
tions.  A  precise  definition  of  all  three  of  these  conditions  as  they 
apply  to  FFT's  is  difficult  to  come  by.  However,  the  effect  of  these 
three  conditions  is  related  to  how  the  FFT  views  the  finite  area  as  one 
period  in  a  domain  infinite  in  extent.  For  a  finite  area  unpadded  with 
zeroes,  there  is  a  chance  that  important  information  clustered  along  the 
edges  can  be  viewed  as  part  of  the  adjacent  period  thus  causing  distor¬ 
tion  in  the  FFT  and  subsequent  image  from  the  IFFT  (11:13,  92-3).  To 
prevent  this  possible  distortion  from  occurring  with  this  project,  the 
original  8x8  array  of  FLIR  data  is  arbitrarily  padded  with  an  addi¬ 
tional  8  rows  of  zeroes  on  the  top,  bottom,  and  sides.  This  additional 
padding  results  in  the  24  x  24  input  array  mentioned  earlier  in  this 
section. 

In  summary,  the  FFT  and  IFFT  processes  for  this  project  is  done 
using  the  FOURT  subroutine.  To  negate  the  translational  effects  on  an 
image,  the  shift  theorem  of  the  Fourier  transform  will  be  digit. illy 
implemented.  Finally,  to  eliminate  possible  distortion  by  edge  effects, 
aliasing  and  leakage  conditions,  the  original  8x8  array  of  data  is 
padded  with  additional  zeroes. 
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In  reality,  the  existence  of  a  noise-free  environment  is  fictitious. 
For  computer  simulations,  this  noise  could  be  added  to  a  computer-gener¬ 
ated  noise-free  image.  The  actual  details  of  how  the  noise  corruption 
simulation  is  accomplished  is  discussed  in  the  next  section.  Nonetheless, 
the  noise  effects  can  be  minimized  by  appropriate  data  processing  tech¬ 
niques.  Since  a  priori  knowledge  of  the  precise  form  of  the  noise  is 
normally  not  available,  the  underlying  mathematics  makes  little  assump¬ 
tion  about  the  precise  form.  However,  it  is  necessary  to  assume  that  the 
noise  changes  faster  from  sample  period  to  sample  period  than  the  image 
pattern  itself.  Traditionally,  the  moving  average  technique  is  used  to 
combat  the  effect  of  noise  under  these  conditions.  This  technique  would 
store  the  most  recent  K-l  pieces  of  data  in  memory  and  average  the  data 
in  memory  with  the  new  data  in  memory  with  the  new  piece  of  data  using  a 
weighting  factor  of  1/K  on  each  piece  of  data  (12:115).  However,  this 
technique  does  have  one  major  disadvantage.  It  would  require  K  storage 

locations  of  computer  memory  for  each  pixel.  Therefore,  for  an  N  x  N 
2 

input  array,  KN  storage  locations  of  computer  memory  is  required.  Thus, 
for  large  K  and  large  input  arrays,  significant  data  processing  problems 
would  arise.  An  alternative  method  which  alleviates  this  problem  associ¬ 
ated  with  the  moving  average  is  called  exponential  smoothing  (12:114-28). 
The  equation  that  is  fundamental  to  this  process  is: 

y (t)  =  a  y ( t)  +  (1-a)  y(t-l)  (12:115)  (16) 

where  y(t)  =  most  recent  averaged  value 
y(t)  =  most  recent  piece  of  data 
9 ( t  —1 )  =  previous  averaged  value 


■ix  =  smoothing  constant 
G  <  (l  <  l 

A  key  parameter  in  the  exponential  smoothing  equation  is  the  smoothing 
constant  a  which  can  vary  anywhere  from  0  to  1  inclusive,  depending  on 
how  much  y(t)  is  to  respond  to  the  most  recent  piece  oi  data.  For  noisy 
data  and  slowly  changing  signal  pattern,  it  is  suggested  that  the  values 
of  a  tend  more  toward  0  than  1  so  that  there  is  some  damping  of  the  noise, 
yet  there  is  still  some  response  to  the  most  recent  piece  of  data  (12: 
115).  For  this  thesis  project,  a  steady  otate  value  of  0.1  for  the 
smoothing  constant  was  assumed.  The  notion  of  steady  state  value  is 
important  since  a  pseudo-exponential  smoothing  technique  is  used  to  pro¬ 
vide  better  sensitivity  to  the  initial  10  pieces  of  data.  In  essence, 
the  smoothing  constant  is  varied  for  each  of  the  first  ten  runs  (i.e. 
a  =  1/K  K  =  1,2, 3... 10)  until  the  steady  state  value  oi  a  =  0.1  was 
reached. 

In  closing  this  section,  it  is  interesting  to  note  that  the  expon¬ 
ential  smoothing  technique  possesses  those  qualities  that  the  moving 

2 

average  lack.  First,  only  2N  computer  storage  locations  are  needed  to 
generate  y(t) .  Second,  y(t)  contains  portions  of  all  past  y(t)'s 
although  the  initial  pieces  of  data  have  less  of  an  impact  on  the  most 
recent  averaged  value. 

Spatial  Noise 

To  model  the  real  world  environment,  noise  corruption  should  be 


added  to  any  computer  simulation.  There  are  two  general  types  of 
characteristics  to  specify  for  the  noise  corruption  in  this  problem, 
temporal  and  spatial.  Temporally  correlated  noise  implies  that  knowing 


something  about  the  noise  at  one  particular  time  infers  something  about 
the  noise  at  a  later  time.  Normally,  this  inference  can  be  modelled  by 
a  correlation  function  that  is  exponential  in  shape  (5:37).  However, 
since  temporally  correlated  noise  corruption  has  little  effect  on  the 
quality  of  the  image  beyond  that  of  temporally  uncorrelated  noise  at  the 
expected  signal-to-noise  ratio  for  this  application,  and  since  it  is 
difficult  to  implement  in  a  computer  simulation,  it  was  not  included  in 
this  thesis  project  (5:37). 

Under  spatial  noise  characteristics,  there  are  two  categories, 
spatially  uncorrelated  and  correlated  noise,  before  examining  these  two 
areas,  a  number  system  for  the  8x8  input  pixel  array  and  the  64  x  64 
i  -sultant  correlation  matrix  have  to  be  developed.  The  8x8  array  is 
numbered  from  1  to  64  by  rows  starting  in  the  upper  left-hand  corner 
(Figure  5) . 
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Figure  5.  1’ixel  Numbering  Scheme  (5:19) 

In  the  associated  64  x  64  correlation  matrix,  each  row  or  column  repre¬ 
sents  how  that  particular  pixel  reLates  to  the  other  pixels  in  the  input 


array.  For  example,  row  1  or  column  1  represents  the  correlation  for 
how  pixel  1  relates  to  itself  and  the  other  63  pixels.  Two  unique 


features  of  the  correlation  matrix  are  (1)  since  the  correlation  between 

two  pixels  is  the  same  in  both  directions  (i.e.  r..  =  r,.),  the  correla- 

ij  Ji 

tion  matrix  is  symmetric  and  (2)  since  the  correlation  coefficient  is 
one  for  a  pixel  related  to  itself,  the  diagonal  elements  of  the  matrix 
are  one.  In  determining  the  off  diagonal  terms,  the  first  and  second 
nearest  neighbor  concept  discussed  in  Harnly  and  Jensen's  thesis  was 
used  (5:18-20).  Essentially,  the  exponential  form  of  the  correlation 
function  is  assumed  in  generating  the  correlation  coefficients.  However, 
the  non  zero  values  generated  are  limited  to  the  first  and  second  nearest 
neighbors  of  the  pixel  in  question.  This  condition  results  in  25  out  of 
64  entires  in  each  row  or  column,  being  nonzero.  Of  these  25  values, 
there  are  only  six  distinct  values;  one  for  the  correlation  coefficient 
of  the  pixel  with  itself,  four  values  of  0.3679,  four  values  of  0.2431, 
four  values  of  0.1353,  eight  values  of  0.1069,  and  four  values  of  0.0591. 
These  nonzero  values,  which  are  used  to  generate  the  64  x  64  correlation 
coefficient  matrix  (£) ,  were  developed  from  the  following  equation 


C.  .  =  exp(—  Ill) 
ij  20 

where  d„  =  distance  from  center  of  pixel  i  to  center  of  pixel  j 


(17) 


The  choice  of  a  correlation  distance  of  20  urad  is  dictated  by  the  non¬ 
zero  values  generated  by  the  first  and  second  nearest  neighbor  concept 
(Figure  6). 
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Figure  6.  First  and  Second  Nearest  Correlations  Coefficient 

For  the  spatially  correlated  case,  the  matrix  of  correlation  coefficients 
becomes 
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To  complete  the  process  for  generating  the  noise,  the  64  x  64  correlation 


coefficient  matrix  is  pro  multiplied  by  the  scalar  value  of  the  variance 
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0^  in  order  to  generate  the  correlation  matrix  iLself. 
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Now  that  _R  is  known,  the  64  x  1  noise  vector  (V^t.)),  which  will  be 
added  to  the  input  array,  can  be  computed  by  using  the  Cholesky  square 
root.  Specifically, 

v/>\  =  9nu/'r\  M  sn 


VitJ  =  7R  W  (t .) 


where  W(t^)  =  64-  dimensional  independent  white  Gaussian  noise  vector 
(zero  mean,  variance  of  I) 

thus  E(U(t.)  WT(t.))  =  1  <5.  . 

-i-i  -  ij 

The  Cholesky  square  root  is  a  unique  matrix  decomposition  which  produces 
the  square  root  of  R  in  the  lower  triangular  form.  This  lower  triangular 
form  minimizes  the  number  of  nonzero  computations  in  generating  V_(t^) . 

To  recover  the  original  _R  matrix,  the  should  be  post-multiplied  by 
fti1  therefore  rl*  -  K.  because  of  this  property,  the  covariance  of 

V(t.)  js  preserved  with  tlio  use  of  the  Cholesky  square  root: 

E(V(t.)  VA(t  ))  =  E(^K  W(t.)  w\t.)  ftl1) 

=  c/l<  Id  ( W  (  C  .  )  w1  ( t  .  )  )  1  =  l/R  1  '/'R1  b.. 

-  -  i  -  j"  —  - -  ij 

=  R  A. .  (19) 

-  J-J 
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As  mentioned  earlier,  the  noise  vector  will  he  added  to  the  noise-free 
input  matrix  to  model  real-world  image  conditions.  This  resulting  vector 
called  z.(t^)  is  the  measurement  vector  that  will  be  provided  to  the 
Extended  Kalman  Filter.  Algebraically,  the  z  vector  is  developed  from 
the  truth  model  and  is  of  the  form 

z(ti)  =  hCxCt^,  ti)  +  V.(ti)  (20) 

where  £(tp  =  output  state  vector  from  the  system  dynamics  model 

h  (x(  t . ) ,  t  J  =  nonlinear  measurement  function  developed  from 
image  intensity  pattern 
V(t.)  =  measurement  noise  vector 

As  discussed  in  the  problem  overview,  the  development  of  the  nonlinear  h 
function  from  the  measurement  vector  z_(t.)  is  a  major  goal  of  this  thesis. 

In  summary,  since  temporally  correlated  noise  has  been  shown  to  have 
little  effect  on  the  Kalman  Filter  performance  at  the  signal-to-noise 
ratios  for  this  project,  only  soatially  correlated  and  uncorrelated  noise 
was  used.  The  correlation  matrix  used  to  generate  the  noise  vector  V(t  ) 
is  developed  from  the  first  and  second  nearest  neighbor  concept  as  it 
applies  to  a  correlation  function  exponential  in  nature.  The  Cholesky 
square  root  of  the  correlation  matrix  is  taken  and  post-multiplied  by  a 
white  noise  vector  W(l^),  with  the  final  result  being  the  noise  vector 
V^(t^)  that  is  added  Lo  the  noise-free  image  array  to  provide  the  neces¬ 
sary  noise  corruption. 

Relationship  between  Truth  Model  and  Extended  Kalman  Filter 

In  the  previous  section,  terms  such  as  truth  model,  system  dynamics, 
and  Extended  Kalman  Filter  were  used.  These  terms  represent  the  heart 
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of  this  projucL.  The  truth  model  is  a  bust  rcprosunlat  ion  ol  the  envir¬ 
onment  from  which  the  measurement  vector  z(t^)  is  generated.  These 
environmental  characteristics  include  Lite  underlying  target  dynamics, 
atmospheric  effects,  and  background  noise.  In  order  to  process  the 
measurement  vector  z_(t^)  and  provide  a  best  estimate  of  the  underlying 
target  centroid  location  for  the  next  sample  period,  the  Extended  Kalman 
Filter  is  used.  The  actual  details  of  the  filter  will  be  discussed  in  a 
later  section.  However,  the  interplay  between  the  Kalman  Filter  and  the 
truth  model  is  a  point  which  at  times  can  be  confusing.  In  order  to 
provide  a  best  estimate  of  the  underlying  centroid  location,  the  dynamic 
models  are  a  facsimile  of  those  in  the  truth  model.  Ideally,  the  models 
should  match  exactly,  but  in  reality,  due  to  the  desire  for  computational 
efficiency,  there  are  always  deviations.  The  amount  of  deviation  in  the 
dynamics  models  for  tile  filter  varies  depending  on  the  problem  applica¬ 
tion.  A  contributing  factor  to  this  deviation  may  be  a  lack  of  knowledge 
concerning  the  exact  form  of  the  truth  models.  However,  due  to  the 
robustness  of  the  Kalman  filter,  adequate  predictions  of  Che  centroid 
location  can  still  be  accomplished. 

Truth  Model 

Much  has  been  discussed  concerning  the  dynamic  models.  For  this 
project,  it  had  been  envisioned  to  investigate  the  feasibility  of  using 
Loth  deterministic  and  stochastic  models.  In  the  realm  of  deterministic 
models,  initially,  the  target  dynamics  involved  a  stationary  target 
centered  in  the  system  FOV.  This  step  was  taken  to  identify  and  elimi¬ 
nate  any  possible  inherent  motion  caused  by  Lhe  pattern  recognition  or 
Extended  Kalman  Filter  algorithms.  The  next  deterministic  model  involved 
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a  target  moving  acorss  the  imago  plane  at  a  constant  velocity.  This  was 
done  to  test  the  robustness  of  the  Extended  Kalman  Filter.  Supposedly, 
the  filter  design  used  for  this  project  can  absorb  any  truth  model  whose 
resultant  motion  from  one  sample  period  to  the  next  is  less  than  one-half 
pixel  (i.e.  less  than  10  urads)  (13). 

In  the  realm  of  stochastic  models,  the  target  and  atmospheric 
dynamics  developed  in  Mercier's  thesis  were  used  (3:9-16).  The  target 
dynamics  are  modelled  in  each  direction  as  a  first-order  Gauss-Markov 


process  driven  by  white  Gaussian  noise 

yt)  =  -J  yon^t)  (16) 

yo  =  -  x  vc)  +  w2(t)  u?) 

where  E(W1(t))  «  E(W2(t))  =  0  (18) 

o 2 

E(W1(t)  W2(s))  =  E(W2(t)  W2(s))  =  3.  6(t-s)  (19) 

AT 

E(W1(t)  W2(s))  =  0  (20) 

X^  =  truth  model  correlation  time 


Vl1(t),  W2(t)  =  continuous  time  independent  white  Gaussian  noise  processes 
2 

=  the  desired  variance  on  the  outputs  X^  and  YQ  (3:10). 

The  atmospheric  jitter  model  was  based  on  a  study  by  the  Analytic 
Sciences  Corporation  (IX; 29, 30)  and  a  data  analysis  by  Hoggc  and  Butts 
(lb).  These  studies  resulted  in  the  development  of  a  third  order  shaping 
filter,  with  a  single  pole  at  14.14  rad/sec  and  a  double  pole  at  659.5 
rad/sec,  that  is  driven  by  white  Gaussian  noise  (3:12). 

Since  any  nC^  order  differential  equation  can  be  written  as  a  set 
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of  coupled  first  order  differential  equations,  the  state  space  model 
incorporating  both  the  target  dynamics  and  atmospheric  effects  was  gener 
ated  using  the  Jordon  canonical  form  (3:13).  Referring  to  the  state 
space  aodel  from  Mercier's  thesis, 

XT  =  F T  XT  +  Gx  Wt  (t)  (21) 

where  £T  =  the  truth  model  plant  matrix 
=  the  truth  model  state  vector 
G.J,  ■  the  truth  model  input  matrix 
W„  =  a  vector  of  white  Gaussian  noise  inputs 

_T 

E (WT ( t) }  =  0  (22) 

E{Wx(t)  WtT(s)}  «  6(t-s)  (23) 

(24) 
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a  =  14.14  rad/sec  b 


659.5  rad/sec 


ab 

(a-b) 


(3:11-14) 


(25) 


(26) 


The  composition  of  these  truth  model  matrices  is  dictated  by  the  compon¬ 
ents  of  the  truth  model  state  vector  describing  the  target  dynamics  .and 
atmospheric  effect  on  the  centroid  location.  For  this  project  from  X^, 

XD  -  X(l)  (27) 


29 


XA  +  X(2)  +  X(3) 
A 


(28) 


=  X(5)  (29) 

Ya  =  X(6)  +  X(7)  (3:14)  (30) 


where  and  are  atmospheric  jitter  variables  as  output  of  third 
order  shaping  filter. 

In  order  to  transition  the  state  vector  from  one  sample  period  to 
another,  propagation  equations  are  developed.  The  basis  for  these 
equations  is  the  solution  to  the  matrix  differential  equation 

±(t,t.)  =  i(t,t.)  (31) 

where  <t>(t,t.)  =  state  transition  matrix  from  time  t.  to  time  t 
—  i  l 

associated  with  the  matrix  JAj. 

♦  (3:14)  (32) 


Since  the  truth  model  plant  matrix  F,p  is  a  constant,  the  state  transi¬ 
tion  matrix  is  stationary  with  respect  to  time.  Therefore,  ^(^..,0 
is  also  constant  for  a  fixed  sample  period  (3:14).  Using  this  truth 
model  plant  matrix,  the  solution  to  the  above  differential  equation 
becomes 


❖x(At) 


-At/XT 
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0 

0 


0 

-a  A  t 

■» 

0 
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0 

-hA.t 


0 

-hAt 

te 

-hAt 


(33) 


At  =  (l 


i+1 


ti> 


The  state  transition  matrix  is  used  in  the  solution  to  the  state  space 
vector  differential  equation 
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St  ‘  4  V°  +  iLf  V°  (34) 

This  solution  is 

VW  =  ^t(Ac)  (35) 

+  £  i+l  iT(ti+1,X)  Gt  (X)  WT(X)dX  (36) 

The  details  for  the  following  argument  are  contained  in  Mercier's  thesis. 
Nonetheless,  the  resulting  equivalent  discrete- time  model  to  the  above 
equation  is 

2Lt  <ti+1>  =  ±T(At)  X^t.)  +  ^Qd  WjCt.)  (37) 

where  is  the  Cholcsky  square  root  of 

Qd  =  /c  i+l  ^(t.^.X)  (^(X)  ^(A)  (^(X)  ^T(ti+1,X)dX  (38) 

B{WJ(ti)}  =  0  (39) 

E{W  (t.)  W,T(t.)}  =  I  A..  (3:14-5)  (40) 

— d  1  — d  1  —  XJ 

In  summary,  three  truth  models  were  developed  for  this  project. 

Two  models  were  deterministic  in  nature;  a  stationary  target  and  a  tar¬ 
get  moving  at  a  constant  velocity  across  the  image  plane.  These  models 
were  primarily  developed  to  trouble-shoot  and  provide  a  bench  mark  for 
the  pattern  recognition  and  Extended  Kalman  Filter  algorithms.  The  third 
model  is  the  stochostic  representation  that  will  ultimately  be  used  to 
analyze  the  Filter  performance. 

Extended  Kalman  Filter 

As  stated  before,  the  basic  Extended  Kalman  Filter  developed  in 
Mercier's  thesis  was  used  for  this  project.  The  target  dynamics  assumed 
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by  the  Filter  is  also  a  first  order  Causs-Markov  process,  generated  by  a 
first-order  lag  driven  by  white  Gaussian  noise  (3:19).  The  differential 
equation  describing  a  particular  state,  which  is  very  similar  to  that  for 
the  stochastic  truth  model,  is 


X(t)  =  -  j  X(t)  +  W(t) 


E(W(t) )  =  0 


(41) 

(42) 


E(W(t)  W(s) )  = 


2  a 


6(t-s) 


(43) 


where  X  ■  correlation  time 
X(t)  =  state 

X(t)  =  time  derivative  of  the  state 
W(t)  =  white  Gaussian  noise  (3:19) 


However,  the  plant,  input,  and  white  noise  covariance  matrices  assumed 

by  the  filter  differ  from  those  developed  in  the  truth  model.  The 

reduced  order  filter  design  assumes  four  states  in  its  state  vector 
T 

(XD,  Yp,  XA>  Yp)  which  results  in  a  state  vector  differential  equation 


X.,(t)  =  Fp  XF(t)  +  WpCt) 


where  X_,(t)  =  filter  state  vector 
r 

F  =>  filter  plant  matrix 
F 

Wp(t)  3  input  white  noise  vector 


7] 


If  = 


(44) 


(45) 


32 


33 


=  i 


(30) 


Since  ^  is  a  time  invariant  matrix,  the  solution  to  the  above  differen¬ 
tial  equation  is  $p(t,t^),  the  filter  state  transition  matrix,  given  by 


V^i)  = 


e^i^D  0 


0 

0 

0 


e-(t'Ci)/AD  0 


0 

0 


-  Ct-t . ) / A 
e  l  A  0 


0 


e  a  A 


(51) 


In  solving  the  associated  differential  equation  for  the  propagation  of 
the  covariance  matrix,  the  following  stochastic  integral  equation  results: 

T, 


-^i+15  “  ^■•(ti+l,Ci)  -(Ci)  (Ci+rti) 


+  /L  i+1  iF(ti+1,X)  iipU)  ^,.T(ci+i,X)ciX 
ci 


(32) 


where  4>„(t .  . ,  ,t .)  =  filter  state  transition  matrix 
— F  l+l  i 

P(t.)  =  conditional  state  covariance  matrix  from  measurement 
—  i 

update  equation  at  time  t^ 

P(t^  )  =  conditional  state  covariance  matrix  propagated  from 

time  t.  to  t.  , 
l  l+l 


Q  (X)  =  noise  covariance  matrix  (3:21) 


Using  the  (^,  matrix  specified  In  equation  ''id,  Lite  solution  to  tin* 
integral  term  in  the  stochastic  equation  above  becomes 
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x(t+)  =  *(t")  +  K(t.)  (z(t.)  -  h(x(t") ,  t.))  (57) 


P(t*)  =  P(t")  -  K(ti)H(ci)P(t")  (58) 

where 


fi 

f! 


P.(t.)  =  propagated  conditional  covariance  matrix  from  filter,  before 
measurement  incorporation  at  time  t^ 

5t(t.)  =  propagated  state  estimate  vector  from  filter 
K(t.)  =  Kalman  Filter  Gain 

—  i 

h(x(t  ),  t.)  =  non-linear  function  of  intensity  measurements  at  time 
“  l  J. 

t.,  as  a  function  of  the  state  estimate 

l 


H(t.) 


3h(x(t-),  t.) 
3x 


linear  function  of  intensity  measurements 
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£(t^)  =  actual  measurement  vector 

To  ease  the  computational  burden  on  the  computer,  a  different  form  of 
the  measurement  update  equation  called  the  inverse  covariance  form  is 
actually  used  in  the  computer  software  (3:26-7).  In  the  conventional 
form,  the  calculation  of  the  Kalman  Filter  Gain  requires  the  inversion 
of  a  64  x  64  matrix  for  each  update,  since  there  are  64  scalar  measure¬ 
ments  (3:26).  However,  in  the  inverse  covariance  form,  only  two  4x4 
inverses  are  performed  to  obtain  the  £(t  )  and  P.(t+)  matrices  (3:26-7). 
The  equations  for  the  inverse  covariance  form  are  as  follows: 


P^U*)  =il_1(t“)  +  HT(t.)K(t.)H(t.)  (59) 

P(t*)  =  (jfV))-1  (60) 

K(t.)  =  P(t^)HT(t.)R"1(t.)  (61) 

«(t*)  =  S(t~)  +  K(t.)  (z(t.)  -  h(x(t.),  t.))  (62) 


Normally,  the  non-linear  h  function  and  its  spatial  derivative 
(linearized  H  function)  which  are  used  in  the  update  equations  are 
explicit  in  form.  However,  for  this  project  the  non-linear  h  function 
is  determined  in  real-time  using  the  FFT,  phase  shifting  and  the  averag¬ 
ing  techniques  discussed  earlier.  While  the  linearized  1!  function  is 
determined  from  the  non-linear  h  function  using  the  numerical  approxima¬ 
tion  discussed  in  Appendix  A. 

In  summary,  the  measurement  update  equation  is  that  portion  of  the 
Kalman  Filter  which  incorporates  the  actual  intensity  measurements  to 
provide  a  best  estimate  of  the  target  centroid  location.  In  order  to 
ease  the  computational  burden  of  inverting  a  64  x  64  matrix,  the  inverse 


covariance 
sion  of  two 


form  of  the  update  equations,  which  only  requires  the  inver- 
4x4  matrices,  was  employed. 
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f.  1 1  Performance  Ann  I  y  s_is 

The  performance  analysts  for  this  thesis  follows  that  outlined  in 
the  Plan  of  Attack  section.  To  begin,  the  validation  of  the  FFT  sub¬ 
routine  FOURT  not  only  demonstrated  its  ability  to  reconstruct  a  two- 

-8 

dimensional  array  of  data,  but  also  established  a  zero  level  of  10 
This  zero  level  represents  the  greatest  non-zero  value  in  the  IFFT  where 
a  zero  would  be  found  in  the  original  image  array. 

The  next  step  involved  the  development  and  validation  of  a  sub¬ 
routine  used  to  perform  the  negating  phase  shift.  This  subroutine  called 
Shift  is  described  in  more  detail  in  Appendix  D.  In  short,  this  sub¬ 
routine  applies  the  complex  conjugate  of  the  resulting  linear  phase  shift 
from  the  translation  of  an  image  in  the  space  domain  to  the  Fourier  trans¬ 
form  of  this  same  image.  As  a  result,  when  the  IFFT  is  taken,  an  image 
which  is  centered  in  the  FOV  is  the  end  product.  Figure  7  illustrates 
this  process  using  the  40  urad  x  40  urad  intensity  pattern. 

In  analyzing  the  performance  of  the  exponential  smoothing  technique, 
a  number  of  parameters  were  varied  using  the  3  gaussian  pattern  as  a  basis 
These  parameters  included  the  signal-to-noise  ratio,  which  varied  between 
10  and  20,  and  the  sigma  values  for  the  spread  of  the  gaussian  patterns 
which  included  1/4  pixel,  1  pixel,  and  3  pixel  spreads.  The  various 
sigma  values  provided  a  spectrum  of  shapes  which  ranged  from  sharply 
peaked  (o  =  1/4  pixel)  to  very  broad  (o  =  3  pixels).  The  actual  results 
are  displayed  in  Figures  8  to  13.  However,  it  is  important  to  note  that 
at  a  signal-to-noise  ratio  of  20,  the  pattern  is  readily  identifiable  ir- 
regardless  of  its  shape.  However,  the  average  technique  does  provide  some 
smoothing  of  the  data.  At  a  signal-to-noise  ratio  of  10,  the  shape  of  the 
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pattern  seems  to  make  a  difference.  For  the  narrow  gaussian  patterns 
(o  =  1/4  pixel)  the  pattern  is  difficult  to  recognize  even  after  50  runs. 
Hut ,  for  the  very  wide  gaussian  pattern  the  general  trend  of  increasing 
and  decreasing  numbers  is  maintained  throughout  the  process,  which  may 
suggest  that  the  filter  performance  may  be  better  for  larger  targets. 
Still,  it  should  be  emphasized  that  this  is  only  a  superficial  analysis 
and  the  more  detailed  analysis  can  only  be  performed  once  a  complete 
filter  simulation  has  been  developed. 

In  implementing  the  first  truth  model,  that  of  a  stationary  target, 

a  major  problem  arose.  Under  ideal  conditions  (i.e.  the  Kalman  filter 

knowing  exactly  where  the  target  is  initially),  the  filter  began  Co 

diverge  from  the  known  target.  This  divergence  initially  appeared  at 

the  first  measurement  update  and  became  progressively  worse  as  each 

subsequent  propagation  and  measurement  update  were  performed.  Initially, 

it  was  speculated  that  there  may  be  a  Kalman  Filter  tuning  problem.  As 

a  result,  various  values  for  the  initial  covariance  matrix  (P_)  to 

—0 

improve  the  initial  response  of  the  filter,  and  for  the  covariance 
matrices  for  the  dynamic  driving  noise  (£)  and  the  measurement  noise  (R) 
to  improve  the  sicady  state  response  were  tried.  Nonetheless,  no  notice¬ 
able  improvement  in  the  filter's  performance  resulted  from  this  investi¬ 
gation.  Next,  it  was  conjectured  that  the  problem  may  be  within  the 
subroutine  UP1)  which  performs  the  Kalman  Filter  measurement  update. 

First,  the  search  for  a  possible  sign  error  or  any  other  fault  in  the 
software  was  performed  without  any  satisfaction.  However,  upon  investi¬ 
gating  tiie  calculation  of  the  Kalman  filter  gain  t<(t.),  it  was  hypothe¬ 
sized  that  the  forward-backwards  difference  method  used  to  calculate  the 
linearized  h  function  (H_(t^))  whicti  in  turn  is  used  to  calculate  the 
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filler  gain,  may  lie  Liu1  sou  ret:  of  the  problem.  Tn  el  imiiiale  any  puss  1 1>  Lo 
concern  that  the  divergence  problem  may  be  caused  by  the  generation  of 
the  unknown  non-linear  and  linear  h  functions,  the  analytic  function  for 
the  noise-free  centered  gaussian  patterns  generated  from  the  subroutine 
Input  3  was  used  as  the  non-linear  h  function.  The  analytic  spatial 
derivative  of  the  gaussian  pattern  in  both  the  x  and  y  directions  was 
next  taken  and  used  as  the  linearized  h  function.  The  results  of  this 
analysis  is  presented  in  Table  1.  As  can  be  seen,  it  is  verified  that 
the  source  of  the  problem  is  not  due  to  the  generation  of  the  unknown  h 
functions.  As  a  possible  solution  to  this  dilemma  involving  the  diver¬ 
gence  of  the  Kalman  Filter,  the  Fourier  transform  derivative  property, 
which  will  be  discussed  in  more  detail  in  the  Conclusion  and  Recommenda¬ 
tion  section,  should  be  implemented. 
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IV  Conclusions  and  Recommendations 


In  conclusion,  this  thesis  did  not  accomplish  nearly  as  much  as 
anticipated  at  the  beginning  of  the  project.  However,  a  number  of  mile¬ 
stones  were  realized.  To  name  a  few,  the  digital  implementation  of  a 
negating  phase  shift,  the  validation  of  the  FOURT  subroutine,  and  the 
implementation  of  the  exponential  smoothing  technique  were  highlights 
of  this  project.  Since  the  investigation  of  the  deterministic  stationary 
truth  model  identified  a  possible  problem  with  the  f orward-backwar is 
difference  approximation,  the  deterministic  constant  velocity  and 
stochastic  truth  models  were  not  analyzed. 

Although  some  analysis  was  accompl ished  with  this  project,  more  is 
needed.  To  begin,  the  arbitrary  assumption  of  padding  the  input  array 
with  8  additional  rows  and  columns  of  zeroes  could  be  inves'. igated  in 
more  detail  to  determine  if  little  distortion  results  in  using  less 
zeroes.  In  addition,  a  further  extension  on  the  analysis  of  the  averag¬ 
ing  technique  would  involve  not  only  a  further  investigation  of  its 
dependency  on  the  image  shape,  but  also  an  investigation  of  how  the 
filter's  performance  is  affected  when  the  steady  state  value  of  the 
smoothing  constant  a  is  varied.  Another  modification  to  the  computer 
algorithm  that  could  be  made  deals  with  the  sampling  scheme  described 
in  Appendix  A.  A  more  unbiased  and  resolvable  sampling  scheme  could  be 
implemented.  However,  the  addition  of  the  Fourier  transform  derivative 
property  would  be  the  most  important  modification  to  the  existing  com¬ 
puter  software,  since  it  provides  the  best  hope  of  solving  the  filter 
divergence  problem.  Analytically,  the  spatial  derivative  of  the  non¬ 
linear  h  function  represents  the  linearized  h  function.  This  spatial 
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derivative  could  be  more  precisely  calculated  by  using  the  following 
Fourier  property: 

T  {— 3^y')-  }  =  +j2TTfx  H(fx,fy)  (63) 

T  {~h-gy>y)  )  -  j 2tt f y  H(fx,fy)  (17:314)  (64) 

As  a  possible  plan  of  attack  for  future  research,  the  following 
steps  may  be  taken.  First  and  most  importantly,  the  filter  divergence 
problem  must  be  solved,  hopefully  with  the  use  of  the  Fourier  Transform 
divergence  property.  Next,  the  deterministic  constant  velocity  and 
stochastic  truth  models  should  be  employed  keeping  in  mind  that  the  nor¬ 
mal  robustness  that  is  characteristic  of  Kalman  filters  may  not  be 
possible  since  the  negating  phase  shift  requires  extremely  accurate 
estimates  of  the  target  location.  Once  this  lias  ben  performed,  the 
filter  performance  can  be  evaluated  based  on  changes  to  the  zero  padding 
of  the  input  array,  changes  to  the  sampling  scheme  of  the  same  input 
array,  and  variation  to  the  steady  state  smoothing  constant. 
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APPENDIX  A 


Numerical  Approxlmat  ion 


la  the  development  of  this  thesis  project,  two  numerical  approxima¬ 
tions  were  used.  In  generating  the  measurement  vector  z^ti),  the  first 
approximation  involved  computing  the  average  intensity  over  each  pixel. 
Algebraically,  the  component  of  the  z(t^)  vector  corresponding  to  a 
particular  pixel  is 


Zjk(ti)  =  1C~  f  ;  ^argot^’V^  +  W 

^  region  of 

..th  .  . 

jk  pixel 


(65) 


where  A  =  area  of  one  pixel 
P 


I(x,y,t^)  =  two-dimensional  intensity  pattern 
v..  (t.)  =  noise  effect  on  the  jkLl’  pixel 

j  k  i 

*  (tt)  =  actual  measurement  from  jktl>  pixel 


Since  a  simple  functional  form  of  tiie  intensity  pattern  that  can  be 
analytically  integrated  is  not  available,  the  exact  integration  of 
equation  65  cannot  be  performed.  Instead,  25  sample  points  from  each 
pixel  will  be  taken  and  averaged  to  provide  the  average  intensity  over 
the  pixel.  Hie  exact  locations  are  given  in  Figure  14  below.  It  should 
be  noted  that  samples  are  taken  along  the  top  and  left  edge  of  the  pixel 
of  interest  so  that  duplicate  sampling  does  not  occur  along  the  edges. 

In  other  words,  the  right  edge  becomes  the  sampled  left  edge  for  the 
next  pixel  to  the  right,  while  the  bottom  edge  becomes  the  sampled  top 
edge  for  the  pixel  just  below.  The  selection  of  this  sampling  scheme 
was  arbitrary.  In  follow-on  projects,  more  samples  from  different  loca¬ 
tions  may  be  taken  to  provide  better  resolution  and  thus  more  sensitivity 
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Figure  14.  Pixel  Sampling  Scheme 

to  the  movement  of  the  intensity  pattern  from  one  sample  period  to 
another. 

The  second  approximation,  which  is  indirectly  tied  to  the  first, 
involves  the  generation  of  the  linearized  h  function.  Since  an  explicit 
analytical  form  of  the  non-linear  h  function  is  not  available,  an  exact 
differentiation  of  this  non-linear  function  to  provide  a  linearized 
function  cannot  be  performed.  Instead,  a  numei’".  ’  approximation  known 
as  the  Forward-Backwards  Difference  Method  is  used  (18:16-17).  A 
familiar  numerical  method  used  in  finite  difference  calculus,  Lhis  tech¬ 
nique  is  used  to  compute  the  derivative  of  a  function  to  any  higher 
order.  However,  to  compute  the  first  order  spatial  derivative  for  a 
particular  pixel,  the  pixel  values,  along  the  direction  for  which  the 
spatial  derivative  is  taken,  just  before  and  jusL  after  the  pixel  of 
interest  is  used.  The  values  are  subtracted  from  each  other  and  divided 
by  the  distance  between  each  other  (i.e.  40  urad) .  For  the  pixels  located 
along  the  edges  for  which  the  pixel  just  before  or  just  afler  does  not 
exist,  either  just  a  forward  or  just  a  backwards  difference  is  taken. 

To  be  more  specific,  the  value  of  the  pixel  itself  is  subtracted  from 

5b 


cither  the  value  of  the  pixel  just  before  or  just  after,  whichever  is 
available,  and  divided  by  the  distance  between  the  two  (i.e.  20  urnd). 
For  example,  referring  to  Figure  5,  Pixel  Numbering  Scheme,  if  the 
spatial  derivative  in  the  x  direction  for  pixel  number  29  is  to  bo  com¬ 
puted,  the  value  for  pixel  28  is  subtracted  from  the  value  for  pixel  30 
and  divided  by  40.  However,  along  the  edge,  just  a  forward  or  backward 
difference  has  to  be  calculated.  For  instance,  to  find  the  spaLisi 
derivative  in  the  y  direction  for  pixel  number  3,  the  value  for  pixel  3 


S 

A 
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is  subtracted  from  the  value  for  pixel  11  and  divided  by  20. 
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Coordinate  System  for  8x8  Input  Array 

The  development  of  a  coordinate  system  to  use  throughout  this  pro¬ 
ject  was  critical.  The  8x8  input  array  represents  an  8  x  8  FLIR  array 
whose  individual  pixel  FOV  is  20  urad  x  20  urad  and  system  FQV  is  160 
urad  x  160  urad.  With  the  additional  8  rows  and  columns  of  zeroes  (see 
Fourier  Transform  section),  the  resulting  24  x  24  array  of  data  actually 
represents  a  coordinate  system  that  is  480  urad  x  480  urad  with  the  (0,0) 
coordinate  in  the  upper  left  hand  corner,  x  increasing  from  left  to  right 
and  y  increasing  from  top  to  bottom  (Figure  15). 


Figure  15.  Coordinate  System  for  Input  Array 
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Input  Patterns 

Although  the  development  of  the  Extended  Kalman  Filter  for  this 
thesis  assumed  no  apriori  knowledge  concerning  the  exact  algebraic  form 
of  the  intensity  pattern,  three  different  intensity  patterns  were  devel¬ 
oped  to  provide  the  measurement  vector  £(t^)  necessary  to  verify  the 
performance  of  the  Kalman  Filter.  The  actual  computer  subroutines 
(Input  1,  Input  2,  Input  3)  are  described  in  Appendix  D.  However,  the 
three  input  patterns  are  a  40  x  40  constant  intensity  block,  3  constant 
height  cylinder  patterns,  and  3  gaussian  profiles.  Figures  16a, b,c 
illustrate  how  the  centered  patterns  would  appear  on  the  FL1R  array. 
Figures  17a, b,c  show  the  resulting  noise-free  values  for  the  average 
intensity  per  pixel.  The  selection  of  these  three  patterns  is  signifi¬ 
cant,  since  the  40  x  40  pattern  provides  a  very  pronounced  difference 
between  zero  and  non-zero  values,  while  the  3  cylinder  pattern  once 
again  provides  a  pronounced  difference  between  zero  and  non-zero  values 
but  is  closer  to  simulating  a  multiple-hot-spot  target,  while  the  3 
gaussian  profiles  provide  the  best  representation  of  the  multiple-hot¬ 


spot  targets. 


Figure  17a.  Discrete  40  x  40  Constant  Intensity 
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Figure  17b.  Discrete  3  Constant  Height  Cylinders 
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Figure  17c.  Discrete  3  Oauss  inn  Profiles 
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